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In  this  report,  we  develop  iterative  algorithms  for  reconstructing 
a  minimum  phase  sequence  from  the  phase  or  magnitude  of  its  Fourier  trans¬ 
form.  The  iterative  techniques  result  in  two  potentially  Important  com¬ 
putational  algorithms.  The  first  is  a  means  of  implementing  the  Hilbert 
transform  of  the  log-magnitude  or  phase  of  the  Fourier  transform  of  a 
minimum  phase  sequence.  This  procedure  avoids  problems  of  phase  tinwrapping 
and  aliasing  Inherent  in  the  direct  discrete  Fourier  transform  implementation 
of  the  Hilbert  transform.  The  second  algorithm  is  a  new  method  of  computing 
samples  of  the  unwrapped  phase.  As  compared  with  other  available  phase 
unwrapping  algorithms,  this  approach  does  not  rely  on  addlng/zifXmultlples 
to  samples  of  the  principal  value  of  the  phase. 
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INTRODUCTION 


Under  certain  conditions  a  signal  can  be  reconstructed  from  a  partial 
specification  in  the  time  domain,  in  the  frequency  domain,  or  in  both 
domains.  A  minimum  or  maximum  phase  signal,  in  particular,  can  be  re¬ 
covered  from  the  phase  or  magnitude  of  its  Fourier  transform  [1].  The 
conventional  reconstruction  algorithm  involves  applying  the  Hilbert  trans¬ 
form  to  the  log-magnitude  or  phase  of  the  Fourier  transform  to  obtain  the 
unknown  component. 

In  this  report,  we  take  an  alternative  approach  by  developing  Iterative 
algorithms  for  reconstructing  a  minimum  (or  maximum)  phase  signal  from  the 
phase  or  magnitude  of  its  Fourier  transform.  Specifically,  we  develop 
algorithms  which  Impose  causality  in  the  time  domain  and  the  given  phase 
or  magnitude  in  the  frequency  domain*  in  an  Iterative  fashion. 

Iterative  algorithms  similar  to  those  we  discuss  here  have  been  useful 
in  a  number  of  areas  where  partial  information  in  the  two  domains  is 
available.  In  particular,  the  algorithms  presented  in  this  paper  are 
similar  in  style  to  the  Gerchberg-Saxton  algorithm  [2],  and  an  iterative 
algorithm  by  Fienup  [3]  in  alternately  Incorporating  partial  information  in 
the  time  and  frequency  domains.  The  Gerchberg-Saxton  algorithm  recovers  a 
two-dimensional  complex  signal  by  iteratively  imposing  the  finite  extent 
of  the  signal  in  the  space-domain  and  its  magnitude  in  both  the  space  and 
frequency  domains.  Similarly,  Fienup' s  algorithm  recovers  a  real  two- 
dimensional  signal  by  iteratively  Imposing  the  finite  extent  and  positivity 
of  the  signal  in  the  space-domain  and  its  magnitude  in  the  frequency  domain. 
Another  iteration  in  this  same  style  recovers  a  finite  length  mixed  phase 
signal  from  the  phase  of  its  Fourier  transform  by  Imposing  a  finite  length 
constraint  in  the  time  domain  and  the  known  phase  in  the  frequency  domain  [4,5]. 
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In  this  report,  we  begin  in  section  2  with  a  discussion  of  a  number 
of  equivalent  conditions  for  a  sequence  to  be  minimum  phase.  In  sections 
III  and  IV,  we  use  these  conditions  in  developing  two  iterative  reconstruc¬ 
tion  algorithms,  one  for  reconstruction  when  the  phase  is  known  and  the 
other  for  reconstruction  when  the  magnitude  is  known. 

In  section  V,  we  discuss  the  discrete  Fourier  transform  (DFT)  reali¬ 
zations  of  the  algorithms  and  illustrate  the  reconstruction  process  with 
examples. 

In  sections  VI  and  VII,  we  describe  two  computational  algorithms 
which  rely  on  the  iterative  procedures  of  sections  III  and  IV.  We  first 
investigate  the  use  of  the  algorithms  in  implementing  the  Hilbert  transform. 
Of  particular  importance  is  reconstruction  of  the  log-magnitude  from  phase 
since  the  iterative  approach  requires  only  the  principal  value  of  the  phase, 
while  the  direct  DFT  implementation  of  the  Hilbert  transform  requires  the 
unwrapped  phase  [6].  The  former  technique,  therefore,  avoids  the  compli¬ 
cations  and  problems  of  phase  unwrapping  which  is  often  computationally 
difficult  [1,7].  In  addition,  the  accuracy  of  the  iterative  process  is  not 
limited  by  the  finite  length  of  the  DFT  as  in  the  direct  approach. 

In  the  second  computational  procedure,  the  Iterative  technique  of 
section  III  is  used  as  the  basis  for  a  new  phase  unwrapping  algorithm. 

This  algorithm  does  not  rely  on  adding  2ir  multiples  to  the  principal  value 
of  the  phase  as  required  by  available  unwrapping  algorithm. 

In  section  VIII,  the  main  results  of  the  report  are  summarized  and 
some  future  research  is  discussed. 
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THE  MINIMUM  PHASE  CONDITION 


In  general,  a  signal  cannot  be  uniquely  specified  by  only  the  phase 
or  magnitude  of  Its  Fourier  transform.  However,  one  condition  under  which 
the  magnitude  and  phase  are  related  is  the  minimum  phase  condition  and 
under  this  condition  a  signal  can  be  uniquely  recovered  from  the  magnitude 
of  its  Fourier  transform  or,  to  within  a  scale  factor,  the  phase  of  its 
Fourier  transform.  In  this  section,  we  discuss  a  number  of  equivalent 
conditions  for  a  signal  to  be  minimum  phase.  These  conditions  will  be  of 
particular  importance  in  section  III  in  developing  the  iterative  algorithms. 

In  the  following  discussion,  we  restrict  the  z-transform  of  the 
sequence  h(n)  to  be  a  rational  function  which  we  express  in  the  form 


(1) 


where  |a^| ,  |b^| ,  jcjJ  and  |d^j  are  less  than  or  equal  to  unity,  zno  is  a 

linear  phase  factor,  and  A  is  a  scale  factor.  When,  in  addition,  h(n)  is 

stable,  i.e.,  E|h(n)|  <  °o,  |c,  |  and  |d.  |  are  strictly  less  than  one. 
n  R  K 

A  complex  function  H(z)  of  a  complex  variable  z  is  defined  to  be 

minimum  phase  if  it  and  its  reciprocal  H  X(z)  are  both  analytic  for 

| z |  1.  A  minimum  phase  sequence  is  then  defined  as  a  sequence  whose 

z-transform  is  minimum  phase.  For  H(z)  rational,  as  in  (2),  the  minimum 

phase  condition  excludes  poles  or  zeros  on  or  outside  the  unit  circle  in 

the  z-plane  or  at  infinity.  As  a  consequence,  the  factors  of  the  form 

(l-bfcz)  corresponding  to  zeros  on  or  outside  the  unit  circle  and  the 

factors  of  the  form  (1-d^z)  corresponding  to  poles  on  or  outside  the  unit 

circle  will  not  be  present.  Furthermore,  in  (2),  n  -0  to  exclude  poles 

o 
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V 


m 


or  zeros  at  infinity.  Thus  for  H(z)  minimum  phase,  (1)  reduces  to 

Mi  -1 

IT  (l-a.  z 

k=l  K 

H(z)  *  A - — 

it1  d-c.z-1) 
k-1  K 

where  |a^|  and  |  CjJ  are  both  strictly  less  than  unity. 

From  (2)  other  conditions  can  be  formulated  for  a  signal  to  be 
minimum  phase.  Two  in  particular  which  we  discuss  below  are  par¬ 
ticularly  useful  in  the  context  Oj.  the  Iterative  algorithms  to  be 
discussed  in  section  III  and  IV. 


Minimum  Phase  Condition  A:  Consider  h(n)  stable  and  H(z)  rational  in  the 
form  of  (1)  with  no  zeros  on  the  unit  circle.  A  necessary  and  sufficient 
condition  for  h(n)  to  be  minimum  phase  is  that  h(n)  be  causal,  i.e., 
h(n)=0  n  <  0,  and  nQ  in  (1)  be  zero. 


From  (2),  it  follows  that  these  conditions  are  necessary.  To  show 
that  they  are  sufficient,  we  want  to  show  that  they  force  (1)  to  reduce  to 
(2).  Clearly,  factors  of  the  form  (1-d^z)  |djJ  <  1  in  the  denominator 
introduce  poles  outside  the  unit  circle  which  would  violate  the  causality 
condition.  With  nQ*0  in  (1),  factors  of  the  form  (1-b^z)  would  introduce 
positive  powers  of  z  in  the  Laurent  expansion  of  H(z),  requiring  h(n)  to 
have  some  non-zero  values  for  negative  values  of  n,  thereby  again  violating 
the  causality  condition.  Therefore,  these  factors  cannot  be  present  and 
with  no«0,  (1)  reduces  to  (2).  Finally,  because  our  condition  assumes  h(n) 
is  stable  and  that  H(z)  has  no  zeros  on  the  unit  circle,  h(n)  is  minimum 
phase . 
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The  above  minimum  phase  conditions  require  that  h(n)  be  causal  and 
that  the  unwrapped  phase  function  have  no  linear  phase  component.  Another 
slightly  different  set  of  necessary  and  sufficient  conditions  for  a  signal 
to  be  minimum  phase  can  be  stated  as  follows: 


Minimum  Phase  Condition  B:  Consider  h(n)  stable  and  H(z)  rational  in  the 
form  of  (1)  with  no  zeros  on  the  unit  circle.  A  necessary  and  sufficient 
condition  for  h(n)  to  be  minimum  phase  is  that  h(n)=0  n  <  0  and  h(0)-A 
where  A  is  the  scale  factor  in  (1). 


Again,  from  (2)  it  follows  that  these  conditions  are  necessary  since 
(2)  has  no  poles  or  zeros  outside  the  unit  circle  or  at  Infinity,  guaranteeing 
causality,  and  from  the  Initial  Value  Theorem  h(0)=lim  H(z)-A.  To  demon- 

Z-*oo 

strate  that  these  conditions  are  sufficient  we  note  that  again  causality 
of  h(n)  will  eliminate  factors  of  the  form  (1-d^z)  In  the  denominator  of 
(1).  Furthermore,  since  the  conditions  require  that  h(n)  be  causal,  the 
Initial  Value  Theorem  can  be  applied  with  the  result  that 


h(0)-  lim  H(z)«lim  A  z  ° 

Z-KXl  Z-KD 


M 

o 

H  (1-b.z) 
k-1  K 


(3) 


Since  h(0)-A, 


lim  z  °  n  (l-b.z)-l  (4) 

z-*»  k-1  k 
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and  since  |bjJ  <1  this  requires  that  nQ=0  and  the  b^'s  be  equal  to  zero. 

Thus,  again  (1)  reduces  to  (2). 

Another  condition  which  can  be  shown  to  be  equivalent  to  minimum  phase 
condition  A  or  B  or  our  original  definition  of  a  minimum  phase  sequence  is 
that  the  log-magnitude  and  unwrapped  phase  of  H(oj)  are  related  through  the 
Hilbert  transform  [1].  The  Hilbert  transform  relation  guarantees  that  a 
minimum  phase  sequence  can  be  uniquely  specified  from  the  Fourier  transform 
phase  and,  to  within  a  scale  factor,  from  the  Fourier  transform  magnitude. 
Consequently,  this  unique  characterization  can  be  made  when  minimum  phase 
condition  A  or  B  holds,  and  the  Fourier  transform  phase  or  magnitude  is  given. 

One  technique  for  minimum  phase  signal  reconstruction  from  phase  or 
magnitude  relies  on  a  DFT  implementation  of  the  Hilbert  transform  [6]. 

Two  drawbacks  to  this  algorithm  are  the  requirement  of  samples  of  the  un¬ 
wrapped  phase  and  inaccuracies  due  to.  aliasing.  In  the  next  two  sections, 
we  develop  iterative  algorithms  for  reconstructing  a  minimum  phase  sequence 
h(n)  from  the  phase  or  magnitude  of  its  Fourier  transform  which  bypass  these 
problems.  Motivated  by  the  minimum  phase  condition  A,  when  the  phase  is 
given  we  impose,  in  an  iterative  fashion,  causality  in  the  time  domain  and 
the  known  phase  (which  has  no  linear  phase  component)  in  the  frequency 
domain.  When  the  resulting  sequence  satisfies  minimum  phase  condition  A 
and  has  the  given  phase,  it  must  equal  h(n)  to  within  a  scale  factor. 

Likewise,  motivated  by  the  minimum  phase  condition  B,  when  the  magnitude 
is  given  we  impose,  in  an  iterative  fashion,  causality  and  the  initial 
value  h(0)  in  the  time  domain,  and  the  known  magnitude  in  the  frequency 
domain.  When  the  algorithm  results  in  a  sequence  which  satisfies  minimum 
phase  condition  B  and  has  the  given  magnitude,  it  must  equal  h(n). 


6 


3 


AN  ITERATIVE  ALGORITHM  FOR  SIGNAL  RECONSTRUCTION  FROM  PHASE 


The  iterative  algorithm  for  reconstructing  a  minimum  phase  signal 
from  its  phase  function  is  shown  in  Figure  1.  The  functions  h^(n),  0^ (gj) 
and  M^(co)  represent  the  signal  and  its  phase  and  magnitude  estimates  on  the 
kth  iteration.  The  function  0^(m)  is  the  known  phase  and  h^^ (n)  -  h^(n)u(n) , 
where  u(n)  is  the  unit  step  function. 

The  algorithm  begins  with  an  initial  guess  M^(u))  of  the  desired  Fourier 
transform  magnitude  and  the  inverse  Fourier  transform  of  (w)  exp  [j  0^(0)) ]  is 

taken.  This  step  yields  h^n),  the  initial  estimate  of  h(n).  Next, 
causality  is  imposed  so  that  h^(n)  is  set  to  zero  for  n<  0  to  obtain 
h^(n).  The  phase  of  the  Fourier  transform  of  h^(n)  is  then  replaced  by  the 
given  phase  and  the  procedure  is  repeated. 

We  now  show  that  the  mean  squared  error  between  h(n)  and  h^(n)  or 
equivalently  between  H(oi)  and  is  non-increasing  on  successive 

iterations.  The  mean  squared  error  on  the  kth  iteration  from  Parseval's 
Theorem  can  be  written  as 


\-irj  l"<“>  -  V“’l: 


=  E  |h(n)  -  h,  (n)  | 


E  |h(n)  -  h,  (n)|  +  E  |h(n) 

n  <  0  n>0 


-  K  (»)| 
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Fig.  1.  Iterative  algorithm  to 
recover  h(n)  from  its  phase. 
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Since  (n)  -  h^n^uCn) , 


|h(n)  -  ^(n)  |  =  |h(n)  -  hk+1(n)|  »  n>°  (6) 


and 


|h(n)  -  hk<n)|  >  |h(n)  -  hfc+1(n)|  =  0  ,  n  <  0  (7) 


Summing  (6)  and  (7)  over  all  n,  we  obtain 

E  «  E  |h(n)  -  h.  (n)  | 2 
k  n  * 

>  E  |hfn)  -  h.  +1(n)|2  (8) 

Next,  from  Parseval's  Theorem,  we  write  (8)  as 

x 

Ek-2F  /  IH(U))  "  \+i(u)f2  dw  (9) 

-tt 

TT 

>  2~~  /|M(W)exp[J0(.)]-Mk+1(.)exp[j0k+1(.)J|2d. 


With  the  triangle  inequality  for  vector  differences,  we  have  at  each 
frequency  w: 

|M(<D)exp[j0(w)]  -  ^+1  (<i))exp[j0k+1(u))]  | 

>  jM(o))exp[j0(w)]  |  -  |Ww)exp[J<Ww)]l 

>  |M(u>)  - 
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Therefore,  from  (9)  and  (10),  and  the  Identity  |exp[j0(u>)] 


2 


Is 


IT 

1  f  2 

Ek  —  2tt  J  |M(u>)  -  \+1^)l  du) 

-7T 

7T  2 

>_  JL_  f  |M(o))exp[j0(a)>]  -  M^CuOexpt  j0(w)]  |  dw 

-7T 

>  ~  J  |H(aj)  -  H|efl0-)|2du) 


-  Ek+1 


(ID 


Since  E^  is,  therefore,  non-increasing  and  has  a  lower  bound  of  zero, 
must  converge  to  a  unique  limit  [8],  The  non-increasing  nature  of  Efc, 
ho-  r,  is  not  sufficient  to  guarantee  that  the  iterates  hfc(n)  converge. 
Neve,  leless,  if  a  converging  solution  with  a  rational  z-transform  exists, 
we  can  show  that 


Aim  h.  (n)  -  ah(n)  (12) 

k-*» 


where  a  is  a  positive  constant 


To  see  this,  note  from  (6),  (7)  and  (10)  that  the  equality  In  (11)  holds 

If  and  only  If  h^(n)  -  h^+1(n)  -  0  for  n<0,  and  8^(o>)  -  0^+1(w). 

Therefore,  since  8^(aj)  contains  no  linear  phase  component  (i.e.,  nQ«0) ,  if 

h^(n)  converges  to  a  sequence  whose  z-transform  is  of  the  form  in  (1),  the 

converging  solution  must  satisfy  the  minimum  phase  condition  A.  Consequently, 

the  converging  solution  is  minimum  phase  with  phase  0.  (w) ,  and  (12)  must  hold.* 

n 

When  h(n)  is  of  finite  duration  (i.e.,  H(z)  has  no  poles),  we  can  impose 
not  only  causality,  but  also  a  finite  duration  constraint  within  the  itera¬ 
tion.  Under  these  particular  constraints,  the  DFT  realization  of  our 
iterative  procedure  (see  section  V)  converges  to  a  limit  of  the  form  in 
(12)  [9]. 


The  constant  a  in  (12)  is  constrained  to  be  positive  since  a  negative  value 
introduces  an  additive  factor  of  tf  into  the  phase  function. 
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AN  ITERATIVE  ALGORITHM  FOR  SIGNAL  RECONSTRUCTION  FROM  MAGNITUDE 


In  this  section  we  present  an  iterative  algorithm  for  reconstruction 
of  a  minimum  phase  signal  from  the  magnitude  of  its  Fourier  transform.  The 
algorithm  is  shown  in  Figure  2.  The  functions  h^n),  6^(01)  and  M^w) 
represent  the  signal,  phase,  and  magnitude  estimates,  respectively  on  the 
kth  iteration  and  h^+^(n)  is  defined  by 


hk+l(n) 


hk(n),  n  >0 
h(0) ,  n  =  0 
0  ,  n  <  0 


(13) 


The  algorithm  begins  with  an  initial  guess  0q(oj)  of  the  desired  phase, 
and  the  inverse  transform  of  M(u))exp[j0p(a)) ]  is  taken  where  M(u>)  «  |H(w)|  is 
the  given  magnitude.  This  step  yields  hp(n),  the  initial  estimate  of  h(n). 
Next,  on  the  basis  of  the  minimum  phase  condition  B,  causality  and  the 
known  value  of  h(0)  are  imposed  so  that  hQ(n)  is  set  to  zero  for  n  <  0 
and  set  to  h(0)  for  n»0,  to  obtain  h. (n).  The  magnitude  of  the  Fourier 
transform  of  h^(n)  is  then  replaced  by  the  given  magnitude  and  the  pro¬ 
cedure  is  repeated. 

It  has  not  been  possible  to  show  that  the  mean  squared  error,  as 
considered  in  section  III,  is  non-increasing  for  this  algorithm.  However, 
an  error  function  that  is  non-increasing  is  the  mean  squared  difference 
between  the  known  magnitude  and  the  estimate  M^w)  on  each  Iteration;  i.e., 

IT 

Ek  “  j  |m(w)  -  M^(to)  |  2du)  (14) 

— TT 
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M(w)exp[j8-(a>)] 


F 

-1 

hk(n) 

f 

Impose  Causality  and 

h(0) 

1 

B 

F 

Vn> 

h(0) 

0 


n  >  C 

n  “ 
n  < 


j  Hfcfl<“)expIjekfl(“)] 

M(W)  \+1(w) 


M(u))exptj9k+1(o))] 


Fig.  2.  Iterative  algorithm  to 
recover  h(n)  from  its  magnitude. 


To  show  that  Efc  is  non-increasing,  we  first  use  the  identity  | exp [ j 9fc (to) ] 1 2-l 
to  express  as 


TT 

\  “  Jn~  J  lM^  _Mk(‘i))|2  |exp[j0k(to)]  |2d(u 

— TT 

TT 

=  2tt~  J  exPf30k(w)]  -  M^ai)  exp fj0k (w)]  |2da) 

-ir 


ir 

“  y*  "  Hk(aj)j2da>  (15) 

-TT 

where  ^(w)  and  H^w)  are  the  Fourier  transforms  of  hfc(n)  and  hfc(n), 
respectively.  From  Parseval’s  Theorem  equation  (15)  is  given  in  the  time 
domain  by 

Ek  =  n  K(n)  "  \(n)*2  (16). 


From  (13),  It  follows  that 


Ih^n)  -  \(n)\  >  ^(n)  -  ^(n)!  -  0  ,  n  >  0 


(17) 


and 

^(n)* «-  h^(fl)  |  -  Ih^n)  -  h^fa)!  *  nl° 
Slimming  (17)  and  (18)  over  all  n,  we  obtain 


(18) 


*  E  Ih^n)  -  hfc(n)|2  >  E  ^(n)  -  h^n)!' 


(19) 


Next,  we  apply  the  triangle  inequality  for  vector  differences,  to  yield 


l\<w)  “  Vi^l  i  "  l\+i(u,)l 


(20) 


Therefore,  we  have  from  Parseval's  Theorem,  and  (19)  and  (20) 


^  >1  Ih^n)  -  \+1(n)|2 


-  -h  J  l"k<")  - 

-TT 


-  17  /*  |M(w)  -  Mk+1(w)|Zdw 

— TT 


!Vi 


(21) 
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Since  is  non-increasing  and  has  a  lower  bound  of  zero,  it  must  converge  to 
a  limit  point  [8]. 

As  with  the  algorithm  in  section  III,  although  we  have  shown  that  the 
error  is  non- increasing ,  we  have  not  shown  that  the  iterates  h^(n)  con¬ 
verge.  However,  if  the  iterates  converge  to  a  sequence  whose  z-transform 
is  rational  with  no  zeros  on  the  unit  circle  and  which  is  causal  with  initial 
value  h(0),  from  the  minimum  phase  condition  B,  the  converging  solution  must 
be  minimum  phase.  Consequently,  if  fn  addition  the  magnitude  of  the  Fourier 
transform  of  the  converging  solution  equals  |h(<d)|,  the  solution  is  the 
unique  minimum  phase  sequence  associated  with  |h(gj)|,  i.e.,  h(n). 

The  convergence  of  h^(n)  has  yet  to  be  rigorously  proven  even  when  a 
finite  length  constraint  is  imposed  within  the  iteration  [9] .  Empirically, 
however,  we  have  found  the  DFT  realization  of  the  algorithm  to  always  con¬ 
verge.  In  the  next  section,  we  shall  illustrate  the  convergence  of  h^(n) 
to  h(n)  with  an  example. 
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5.  REALIZATIONS  OF  THE  ITERATIVE  ALGORITHMS  USING  THE  DFT 

Since  the  iterative  algorithms  will  be  implemented  on  a  digital  computer, 
we  can  compute  a  Fourier  transform  at  only  a  finite  number  of  points.  In 
particular,  we  shall  use  the  DFT. 

One  consequence  of  the  DFT  realization  is  that  our  desired  sequence  h(n) 
must  be  of  finite  duration.  Imposing  a  finite  duration  constraint  within  the 
iterations,  however,  does  not  change  the  non- increasing  nature  of  the  error 
functions,  as  can  be  seen  from  (8)  and  (19). 

A  second  consequence  of  the  DFT  realization  is  that  only  uniformly 
spaced  samples  of  the  phase  and  magnitude  functions  are  available.  Neverthe¬ 
less,  it  is  again  possible  to  show  that  the  non-increasing  nature  of  E^  is 
not  altered  when  we  use  samples  of  the  magnitudes  and  Fourier  transforms  in 
the  expressions  for  E^  in  (5)  and  (14)  [4,10]. 

Finally,  questions  of  convergence  need  to  be  addressed.  Consider  first 
minimum  phase  reconstruction  from  phase  samples.  When  H(z)  is  constrained 
to  have  no  conjugate  reciprocal  zero  pairs  and  no  zeros  on  the  unit  circle, 
a  unique  sequence  h(n)  of  length  M  (to  within  a  scale  factor)  is  guaranteed 
when  given  M-l  or  more  phase  samples  of  0^(w)  in  the  open  frequency  interval 
(0,ir)  [4].  A  minimum  phase  sequence,  in  particular,  satisfies  these  con¬ 
straints.  Therefore,  the  DFT  realization  of  the  iterative  algorithm  to 
reconstruct  a  minimum  phase  sequence  from  its  phase  samples  can  be  implemented 
with  a  DFT  of  length  N^2M.  Furthermore,  this  iteration  will  converge  to 
ah(n)  for  0<^  n  £N-1,  where  a  is  positive  [9]. 

Consider  next  the  dual  problem  of  developing  a  DFT  realization  of  the 
iterative  algorithm  to  recover  a  minimum  phase  sequence  of  length  M  from  a 
magnitude  function.  In  this  case,  there  exists  only  one  M-polnt  sequence, 
l.e.,  the  minimum  phase  sequence  h(n),  when  h(0)  is  specified  along  with  M 
or  more  uniformly  spaced  samples  of  the  magnitude  in  the  half  open  frequency 
interval  f 0,tt)  [10].  Therefore,  a  DFT  realization  of  the  iteration  can  be 
implemented  with  DFT  length  N  >_  2M-1.  If  the  algorithm  converges  to  a  causal 
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sequence  of  length  M  with  initial  value  h(0)  and  the  known  magnitude  samples, 
the  converging  solution  must  equal  h(n)  for  0£n  £N-1. 

We  now  consider  two  examples  where  the  DFT  length  is  512  points  which 
is  twice  the  length  of  h(n).  In  the  first  example,  the  initial  magnitude 
guess  is  unity,  and  in  the  second  example  the  Initial  phase  guess  is  zero. 

Example  5.1:  Signal  Reconstruction  From  Phase 

Consider  a  256-point  minimum  phase  signal  h(n)  illustrated  in  Figure  3. 
The  phase  is  known  and  we  wish  to  reconstruct  h(n).  The  functions  h^(n) 
and  log[M^(ui)]  are  depicted  in  Figures  3  and  4  along  with  the  originals  for 
k  equal  to  1,  5  and  45.  The  signal  h^(n)  (to  within  a  multiplicative  constant) 
and  the  spectrum  loglM^Cw)]  (to  within  an  additive  constant)  are  indistintuish- 
able  from  the  originals  after  45  iterations. 


Example  5.2:  Signal  Reconstruction  From  Magnitude 

In  this  example,  we  consider  the  sequence  of  example  5.1,  but  where 
the  Fourier  transform  magnitude  is  given.  The  functions  h^(n)  and  0^(0)) 
are  depicted  in  Figure  5  and  6  with  the  originals  for  k  equal  to  1,  5  and 
25.  The  functions  h^(n)  and  0^(0))  are  indistinguishable  from  the  originals 
after  25  iterations. 


(c)  5  iterations 


(d)  25  iterations 


Fig.  6.  Continued 


6.  IMPLEMENTATION  OF  THE  HILBERT  TRANSFORM 


In  this  section  and  the  next,  we  Investigate  two  computational  algorithms 
based  on  the  iterative  procedures  of  sections  III  and  IV:  (i)  a  new  means 
of  implementing  the  Hilbert  transform,  and  (ii)  the  use  of  this  implementation 
as  the  basis  for  a  new  phase  unwrapping  algorithm. 

For  a  minimum  phase  signal,  the  log-magnitude  and  phase  of  the  Fourier 
transform  are  related  through  the  Hilbert  transform  and  the  direct  imple¬ 
mentation  of  the  Hilbert  transform  using  the  DFT  has  been  extensively  investi¬ 
gated  [6].  One  disadvantage  of  this  implementation  is  that  in  computing  the 
log-magnitude  from  the  phase,  samples  of  the  unwrapped  phase  are  required 
and  are  often  difficult  to  compute.  A  second  drawback  is  that  aliasing 
occurs  in  the  inverse  discrete  Fourier  transform  of  samples  of  the  log- 
magnitude  and  unwrapped  phase  due  to  a  finite  DFT  length,  limiting  the 
accuracy  of  the  computed  samples  of  the  unknown  component. 

An  alternative  to  the  direct  implementation  of  the  Hilbert  transform 
exploits  the  iterative  algorithms  of  sections  III  and  IV.  When  the  phase 
is  given,  through  the  use  of  the  algorithm  in  section  III,  ah(n)  is  first 
obtained  from  the  phase  and,  in  particular,  does  not  require  samples  of  the 
unwrapped  phase.  From  ah(n)  the  log-magnitude  of  ah(co),  representing  the 
Hilbert  transform  of  the  phase  to  within  an  additive  factor  is  then  computed. 
Furthermore,  by  increasing  the  number  of  iterations  we  can  come  arbitrarily 
close  to  samples  of  the  log-magnitude  so  that  accuracy  is  not  limited  by  a  fixed 
DFT  length  as  in  the  direct  approach. 

A  similar  procedure  can,  of  course,  be  applied  through  the  use  of  the 
iterative  algorithm  in  section  IV  to  implement  the  Hilbert  transform  of  a 
given  logmagnitude  function.  As  before,  for  a  fixed  DFT  length,  we  can  come 
arbitrarily  close  to  samples  of  the  phase  by  increasing  the  number  of  itera¬ 
tions.  If  h(0)  is  not  known  a  priori  (recall  (13)),  it  can  be  obtained 
(at  least  in  theory)  from  the  magnitude,  although  in  practice  h(0)  can  be 
computed  only  approximately  [1],  However,  it  was  found  empirically  that  the 
iterates  always  converge  to  h(n)  when  only  causality  is  imposed  in  the  time 
domain  (i.e.,  h(0)  is  assumed  unknown)  and  the  initial  phase  0q (go)  is  set  to  zero. 
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7.  A  NEW  PHASE  UNWRAPPING  ALGORITHM 


There  are  a  variety  of  applications  in  which  it  is  desired  to  obtain  an 
unwrapped  phase.  Current  algorithms  rely  on  adding  multiples  of  2tt  to 
samples  of  the  principal  value  of  the  phase  [1,7].  In  this  section,  we 
present  a  phase  unwrapping  algorithm  which  avoids  such  considerations  and 
which  relies  primarily  on  the  iterative  algorithm  of  section  III. 

We  assume  either  that  there  is  no  linear  phase  component  in  H(z)  as 

n 

given  in  (2),  or  that  we  can  separately  determine  z  °.  Let  0(io)  denote  the 
desired  unwrapped  phase  function  and  6  (to)  its  value  modulo  2ir.  The  proposed 
phase  unwrapping  algorithm  proceeds  as  follows: 


(i)  Remove  the  linear  phase  component  to  obtain 

the  principal  value  of  the  phase  of  H(u>)exp[-jn  to], 
denoted  by  ©^((o). 

(ii)  Apply  the  iterative  algorithm  of  section  III 

ylth  a  causality  constraint  and  with  phase 

0  (to)  to  obtain  a  minimum  phase  sequence  h  (n) . 

P  mp 

(iii)  Compute  loglH  (co)  I  where  H  (to)  is  the  Fourier 

mp  mp 

transform  of  h  (n) , 
mp 

(iv)  Apply  the  Hilbert  transform  to  log|H  (co)  |  to 
obtain  the  unwrapped  phase  function  0(<o)-noio. 

(v)  Add  back  the  linear  phase  component  to  obtain 
the  desired  unwrapped  phase. 


Of  particular  interest  is  step  (ii)  which  yields  the  same  minimum  phase 

sequence  h  . (n)  that  would  be  obtained  by  a  Hilbert  transform  of  the  un- 
mp 

wrapped  phase  but  bypasses  the  need  for  phase  unwrapping. 


There  are  two  major  considerations  in  the  use  of  this  algorithm.  First, 

the  minimum  phase  sequence  h  (n)  derived  from  the  iteration  is  of  infinite 

mp 

extent  regardless  of  whether  the  original  sequence  h(n)  is  of  finite  duration 

[10].  Therefore,  a  possible  problem  with  aliasing  arises.  The  DFT  length 

must  be  sufficiently  large  so  that  the  minimum  phase  sequence  h  (n)  decays 

mp 

effectively  to  zero.  In  particular,  when  h  (n)«0  for  n  >M  the  DFT  length, 

mp 

from  our  discussion  in  section  V,  should  be  at  least  2M. 

The  second  consideration  is  the  linear  phase  factor  of  H(z) .  The 
presence  of  this  term  represents  a  potential  drawback  to  the  algorithm  since 
a  priori  knowledge  of  such  a  factor  is  sometimes  difficult  to  obtain. 


8.  SUMMARY  AND  CONCLUSIONS 

In  this  report,  we  have  developed  iterative  algorithms  for  recon¬ 
structing  a  minimum  phase  sequence  from  either  the  phase  or  magnitude  of 
its  Fourier  transform.  When  the  phase  is  known,  the  mean  squared  error 
between  the  desired  Fourier  transform  and  its  estimate  was  shown  to  be  non¬ 
increasing  on  successive  iterations.  Likewise,  when  the  magnitude  is 
given,  on  successive  iterations,  the  mean  squared  error  between  the  known 
magnitude  and  its  estimate  is  non-increasing.  In  addition,  we  noted  that 
convergence  of  the  iteration  with  known  phase  samples  (i.e.,  the  DFT 
realization)  has  been  demonstrated,  but  convergence  of  the  iteration  with 
magnitude  samples  has  been  shown  only  empirically. 

Finally,  we  developed  two  computational  algorithms  based  on  the 
iterative  procedures:  (i)  a  new  means  of  implementing  the  Hilbert  transform 
which  doesn't  require  an  unwrapped  phase  and  potentially  provides  greater 
accuracy  than  the  direct  approach,  and  (ii)  a  new  phase  unwrapping 
algorithm  which  doesn't  require  adding  multiples  of  2tt  to  the  principal 
value  of  the  phase. 

The  iterative  algorithms,  as  presented,  rely  on  exact  knowledge  of 
the  magnitude,  phase,  and  the  initial  value  of  the  desired  signal.  Sensi¬ 
tivity  to  the  inexactness  of  these  quantities,  to  quantization  noise,  and 
other  forms  of  degradation  is  not  understood  and  is  a  significant  area 
of  future  research. 

In  practice,  we  have  found  that  the  iterative  algorithm  converge 
sometimes  slowly  (e.g.,  after  several  hundred  iterations)  and  sometimes 
quickly  (e.g.,  after  a  few  iterations).  Consequently,  determining  rates  of 
convergence  in  terms  of  characteristics  of  the  minimum  phase  signal  and 
initial  magnitude  or  phase  estimates,  and  methods  of  speeding  up  convergence 
should  be  explored. 

Another  area  being  considered  is  the  interchange  of  the  signal  re¬ 
construction  problems.  In  particular,  we  have  found  empirically  that 
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when  |h((u)|  and  6^(u)  are  effectively  interchanged  through  j£ogH(u),  a 
slightly  modified  version  of  the  iterative  algorithm  of  section  III, 
requiring  a  phase  function,  will  recover  h(n)  from  its  magnitude.  Likewise, 
when  the  phase  is  known,  h(n)  is  recovered  by  a  procedure  similar  to  the 
Iteration  in  section  IV  which  requires  a  magnitude  function.  These  results 
have  led  to  some  interesting  theoretical  speculations  about  the  duality 
of  the  reconstruction  problems  and  their  iterative  solutions. 


REFERENCES 


(1)  A.V.  Oppenheim  and  R.W.  Schafer,  Digital  Signal  Processing  (Prentice- 
Hall,  Englewood  Cliffs,  New  Jersey,  1975). 

(2)  R.W.  Gerchberg  and  W.O.  Saxton,  "A  Practical  Algorithm  for  the 
Determination  of  Phase  From  Image  and  Diffraction  Plane  Pictures," 
Optik  35,  237  (1972). 

(3)  J.R.  Fienup,  "Reconstruction  of  an  Object  from  the  Modulus  of  its 
Fourier  Transform,"  Opt.  Lett.  3^,  27  (1978). 

(4)  M.  Hayes,  J.S.  Lim,  and  A.V.  Oppenheim,  "Signal  Reconstruction  from 
Phase  and  Magnitude,"  submitted  to  IEEE  Trans.  Acoust.,  Speech, 
and  Signal  Processing. 

(5)  A.V.  Oppenheim,  M.H.  Hayes,  and  J.S.  Lim,  "Iterative  Procedures  for 
Signal  Reconstruction  from  Phase,"  SPIE  Proceedings,  April  1980 

(to  be  published) . 

(6)  B.  Gold,  A.V.  Oppenheim,  and  C.M.  Rader,  "Theory  and  Implementation 
of  the  Discrete  Hilbert  Transform,"  Proceedings  of  the  Symposium  on 
Computer  Processing  in  Communications,  Vo 1.  19  (Polytechnic  Press, 

New  York,  1970). 

(7)  J.M.  Tribolet,  "A  New  Phase  Unwrapping  Algorithm,"  IEEE  Trans. 

Acoust.,  Speech,  and  Signal  Processing  ASSP-25.  170  (1977). 

(8)  M.  Rosenlight,  Introduction  to  Analysis  (Scott,  Foresman,  Glenview, 
Illinois,  1970). 

(9)  V.T.  Tom,  T.F.  Quatieri,  M.H.  Hayes,  and  J.H.  McClellan,  "Convergence 
of  Iterative  Nonexpansive  Signal  Reconstruction  Algorithms,"  Technical 
Note  1980-21,  Lincoln  Laboratory,  M.I.T.  (1  August  1980). 

(10)  T.F.  Quatieri,  "Phase  Estimation  with  Application  to  Speech  Analysis  - 
Synthesis,"  Sc.D.  Thesis,  Dept.  Elec.  Eng.  and  Computer  Science,  MIT, 
Cambridge,  MA  (November  1979). 


32 


'Iterative  Techniques  for  Minimum  Phase  Signal 
Reconstruction  from  Phase  or  Magnitude  p 


trmrommtt  555nilWi»T  nSmSBT 
Technical  Note  1980-34 


7.  AUTHORS 


I.  CONTRACT  OR  GRANT  NUMBERTU 


A)  Thomas  FyQuatieri^^ Alan  V^Oppenheim  yt 


5 


NAG5-2 

— FI952S-83-C-0K|O2' 
N00014-75-C-JW5 


+- . -i-a _ 


R049-328 


».  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Lincoln  Laboratory,  M.I.T.  M.I.T. 

Research  Laboratory  of  Electronics, 

n _ 1 _ u. _ w  i  Anson  '  / 


P.O.  Box  73 
Lexington,  MA  02173 


Cambridge,  MA  02139 


10.  PROGRAM  ELEMENT,  PROJECT,  TASK 
AREA  A  WORK  UNIT  NUMBERS 

Program  Element  No .  62702  F 
Project  Nor^.  4594  / 


11.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Air  Feeee  System  Command,  U5AF 
•  An 

k  DC  20331 


Goddard  Space  FU#t  Conner 
Cedi  932 

GreetfMlt.  MD  20771 


Hi 


August  198)/  "V 


to  MiiyirP  ftp  piftH 
••s  imm  i  Rvtt 

40 


14.  MONITORING  AGENCY  N  AME  A  ADDRESS  (if  different  from  CoraroUin,  Office) 


IS.  SECURITY  CLASS,  (of  tkUttport) 

SdA 


Electronic  Systems  Division 
Hanscom  AFB 
Bedford,  MA  01731 


Office  of  Naval  Research 
800  N.  Quincy  Street 
Arlington,  VA  22217 


Unclassified 


ISa.  DECLASSIFICATION  DOWNGRADING 
SCHEDULE 


M.  DISTRIBUTION  STATEMENT  (of  tbit  Report) 


Approved  for  public  release;  distribution  unlimited. 


17.  DISTRIBUTION  STATEMENT  (of  Ido  (Attract  enured  <11  Block  ]0,  If  different  pom  Report ) 


IB.  SUPPLEMENTARY  NOTES 


None 


It.  KEY  WORDS  (Continue  on  reverse  tide  if  n 

signal  reconstruction 

phase 

magnitude 


eery  end  identify  by  block  number) 

minimum  phase  sequence 
iterative  algorithms 
Fourier  transform 


Hilbert  transform 
time  domain 
frequency  domain 


SO.  ABSTRACT  (Continue  a*  reverie  tide  if  neeeeetry  end  identify  by  Mock  number) 

In  this  report,  we  develop  Iterative  algorithms  for  reconstructing  a  minimum  phase  sequence  from 
the  phase  or  magnitude  of  its  Fourier  transform.  The  iterative  techniques  result  in  two  potentially  im¬ 
portant  computational  algorithms .  The  first  is  a  means  of  Implementing  the  Hilbert  transform  of  the 
log-magnitude  or  phase  of  the  Fourier  transform  of  a  minimum  phase  sequence.  This  procedure  avoids 
problems  of  phase  unwrapping  and  aliasing  Inherent  in  the  direct  discrete  Fourier  transform  implemen¬ 
tation  of  the  Hilbert  transform .  The  second  algorithm  is  a  new  method  of  computing  samples  of  foe 
unwrapped  phase.  As  compared  with  other  available  phase  unwrapping  algorithms,  this  approach  does 
not  rely  on  adding  2*  multiples  to  samples  of  the  principal  value  of  foe  phase. 


_ UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  OF  THIS  FA0C  (When  thee  toured) 

ion  so 

- - - — - - — 


DO  ,  J^7,  1471  EWTION  OF  1  NOV  «S  IS  OBSOLETE 


